#---- Replication script 4/4  ----

#Turnbull-Dugarte
#"Liberal in the sheets [...]"

###This script produces descriptive analysis from YouGov data and the meta-comparison of gender on LGBTQ+ attitudes from the ESS

remove(list=ls())

#install.packages("pacman")

library(pacman)
p_load(here, tidyverse, dplyr, brooom, haven, modelsummary, rcartcolor, jtools,
       list, ggpubr, sruvey, xtable, interactions, viridisLite,
       ggdist, patchwork, estimatr, margins, grf, starbility, miceadds)


###FIGURE 1###
###BES data not included in rep. files due to dataing sharing regulations but accessible from: 
### https://www.britishelectionstudy.com/data/
col3<- c("#D60270", "#9b4f96", "#0038a8")

df_BES <- read_dta("data/BES_W25.dta")
df_BES<- df_BES  %>% 
  mutate(age_factor=as.factor(age_factor),
         sexuality_factor=as.factor(sexuality_factor),
         gender_factor=as.factor(gender_factor))

df_prop <- df_BES %>%
  group_by(age_factor, sexuality_factor) %>%
  summarise(count = n()) %>%
  mutate(proportion = count / sum(count)) %>%
  ungroup()
df_prop <- df_prop %>% 
  filter(sexuality_factor!="Straight")

df_prop <- df_prop %>%
  mutate(label = paste0(round(proportion * 100, 1), "%"))

plot_full <- ggplot(df_prop, aes(x = age_factor, y = proportion, fill = sexuality_factor)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
  geom_label(aes(label = label, group = sexuality_factor, color = sexuality_factor),
             position = position_dodge(width = 0.9),
             hjust = -0.1, size = 2.5, fill = "white") +
  scale_fill_manual(values = col3) +
  scale_color_manual(values = col3) +
  labs(title = "Sexual identities across cohorts",
       x = "", y = "", fill = "Sexuality") +
  theme_ggdist() +
  scale_y_continuous(
    limits = c(0, 0.2),
    labels = scales::percent_format(scale = 100)) +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold")) +
  annotate("label", x = 3, y = .1, 
           label = "Bisexual", size = 3, color = "#D60270", fontface = "bold", hjust = 0) +
  annotate("label", x = 3.4, y = .1, 
           label = "Gay/Lesbian", size = 3, color = "#9b4f96", fontface = "bold", hjust = 0) +
  annotate("label", x = 3.8, y = .1, 
           label = "Other", size = 3, color = "#0038a8", fontface = "bold", hjust = 0) +
  coord_flip()

df_prop_gender <- df_BES %>%
  group_by(gender_factor, age_factor, sexuality_factor) %>%
  summarise(count = n()) %>%
  mutate(proportion = count / sum(count)) %>%
  ungroup()
df_prop_gender <- df_prop_gender %>% 
  filter(sexuality_factor!="Straight")
df_prop_gender <- df_prop_gender %>%
  mutate(label = paste0(round(proportion * 100, 1), "%"))

plot_gender <- ggplot(df_prop_gender, aes(x = age_factor, y = proportion, fill = sexuality_factor)) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9)) +
  geom_label(aes(label = label, group = sexuality_factor, color = sexuality_factor),
             position = position_dodge(width = 0.9),
             hjust = -0.1, size = 2.5, fill = "white") +
  scale_fill_manual(values = col3) +
  scale_color_manual(values = col3) +
  labs(title = "Gender differences in sexual identities across cohorts",
       x = "", y = "", fill = "Sexuality",
       caption = "Data: British Election Study (N:30,363)") +
  scale_y_continuous(
    limits = c(0, 0.3),
    labels = scales::percent_format(scale = 100)) +
  theme_ggdist() +
  theme(legend.position = "none",
        strip.background = element_blank(),
        plot.title = element_text(face = "bold"),
        strip.text = element_text(face = "bold", size = 12, hjust=0)) +
  facet_wrap(~ gender_factor) +
  coord_flip()

plot_full/plot_gender
ggsave("figures/Figure1.png",  height=7, width=7, dpi=600)
ggsave("figures/tifs/Figure1.tif",  height=7, width=7, dpi=600)

###FIGURE 2###
###ESS data not included in rep. files due to dataing sharing regulations but accessible from: 
### https://ess.sikt.no/en/
df <- read_dta("data/ESS1_10C.dta")
df <-as.data.frame(df)
df<- df  %>% 
  mutate(countryyear=as.factor(countryyear),
         cntry=as.factor(cntry),
         essround=as.factor(essround))

df2 <- df %>%
  filter(!is.na(hmsacld_std_reverse))
df3 <- df %>%
  filter(!is.na(hmsfmlsh_std_reverse))

fit_model1 <- function(data) {
  if (is.data.frame(data)) {
    data <- droplevels(data)  
  }
  lm_robust(freehms_std_reverse ~ gndr, data = data)
}

fit_model2 <- function(data) {
  if (is.data.frame(data)) {
    data <- droplevels(data)  
  }
  lm_robust(hmsacld_std_reverse ~ gndr, data = data)
}

fit_model3 <- function(data) {
  if (is.data.frame(data)) {
    data <- droplevels(data)  
  }
  lm_robust(hmsfmlsh_std ~ gndr, data = data)
}



model_results1 <- df %>%
  group_by(countryyear) %>%
  do(tidy(fit_model1(.)))

model_results2 <- df2 %>%
  group_by(countryyear) %>%
  do(tidy(fit_model2(.)))

model_results3 <- df3 %>%
  group_by(countryyear) %>%
  do(tidy(fit_model3(.)))


coefficients1 <- model_results1 %>%
  filter(term == "gndr") %>%
  select(countryyear, estimate, conf.low, conf.high)
coefficients1 <- coefficients1 %>%
  mutate(outcome = "one")

coefficients2 <- model_results2 %>%
  filter(term == "gndr") %>%
  select(countryyear, estimate, conf.low, conf.high)
coefficients2 <- coefficients2 %>%
  mutate(outcome = "two")

coefficients3 <- model_results3 %>%
  filter(term == "gndr") %>%
  select(countryyear, estimate, conf.low, conf.high)
coefficients3 <- coefficients3 %>%
  mutate(outcome = "three")

combined_coefs <- bind_rows(coefficients1, coefficients2, coefficients3)


combined_coefs <- combined_coefs %>%
  mutate(coef = paste(countryyear, outcome, sep = "_"))
combined_coefs <- combined_coefs %>%
  mutate(sig = ifelse(estimate > 0, "yes", "no"))
table(combined_coefs$sig)
colors <- c("grey70","black")


plot<- ggplot(combined_coefs, aes(x= reorder(coef, estimate), y = estimate, color=sig)) +
  geom_hline(yintercept = 0, linetype="longdash", colour="red3")+
  scale_color_manual(values = colors, guide = "none")+
   geom_point() +  
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high)) +  
  labs(x = "Individual outcome-specific country-year coefficient", 
       y = "Bivariate difference in means (women vs. men)",
       title="Gender gap in LGBTQ+ tolerance across Europe",
       subtitle="Observations (N=425): gender comparisons across LGBTQ+ response items, countries, & years") +  
  theme_ggdist()+
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        plot.title=element_text(color="black", face="bold"),
        legend.position = "none")
ggsave("figures/Figure2.png", dpi=600, height=5, width=8)
ggsave("figures/tifs/Figure2.tif", height=5, width=8, dpi=600)


###FIGURE D1###
df_ipsos <- read_dta("data/misc/ipsos_props.dta")
df_ipsos<- df_ipsos  %>% 
  mutate(country=as.factor(country),
         highlight=as.factor(highlight))
highlight<- c("grey90", "#C5407EFF")
highlight<- c("#C5407EFF")


highlight<- c("grey90", "#C5407EFF")
ggplot(df_ipsos, aes(x = reorder(country, prop), y = prop, fill = highlight)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values=highlight)+
  labs(title = "Net percent (%) of LGBTQ identifying adults",
       caption = "Data: Country-level estimates from Ipsos Pride Report (2023)",
       subtitle = "Thirty-country average:9%",
       x = "",
       y = "")+
  theme_ggdist()+
  ylim(0,.17)+
  theme(legend.position = "none",
        plot.title = element_text(face="bold"),
        axis.text.y = element_text(face="bold"))+
  geom_hline(yintercept=.09, linetype="dashed")+
  scale_y_continuous(
    labels = scales::percent_format(scale = 100))+
  coord_flip()
ggsave("figures/appendix/FigureD1.png", dpi=500)


###FIGURE D2###
col3<- c("#D60270", "#9b4f96", "#0038a8")
df_yougov <- read_excel("data/misc/yougov_3way.xlsx")
df_yougov_long <- df_yougov %>%
  pivot_longer(cols = -age, names_to = "sexuality", values_to = "value")
df_yougov_long <- df_yougov_long %>% 
  filter(sexuality!="Straight")


plot_yougov<- ggplot(df_yougov_long, aes(x = age, y = value, fill = sexuality)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values=col3)+
  labs(title = "Non-heterosexual sexual identities across cohorts",
       caption = "Data: Yougov UK tracker (August 2023)",
       x = "",
       y = "",
       fill = "Sexuality") +
  theme_ggdist()+
  ylim(0,.5)+
  scale_y_continuous(
    labels = scales::percent_format(scale = 100))+
  theme(legend.position = "none",
        plot.title = element_text(face="bold"))+
  annotate("label", x = 3.2, y = .3, 
           label = "Asexual", size=4, color="#D60270", fontface="bold", hjust=0)+
  annotate("label", x = 3.4, y = .3, 
           label = "Between straight & Gay/Lesbian", size=4, color="#9b4f96", fontface="bold", hjust=0)+
  annotate("label", x =3.6 , y = .3, 
           label = "Gay/Lesbian", size=4, color="#0038a8", fontface="bold", hjust=0)+
  coord_flip()
ggsave("figures/appendix/FigureD2.png",dpi=500)


###FIGURE D3###
df_yougov2 <- read_excel("data/misc/yougov_scale.xlsx")
df_yougov_long2 <- df_yougov2 %>%
  pivot_longer(cols = -gender, names_to = "view", values_to = "value")
df_yougov_long2<- df_yougov_long2  %>% 
  mutate(value=as.numeric(value))

plot_yougov_scale<- ggplot(df_yougov_long2, aes(x = view, y = value, fill = gender)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values=col3)+
  labs(title = "Thinking about sexuality, which of the following comes closer to your view?",
       caption = "Data: Yougov UK tracker (August 2023)",
       x = "",
       subtitle="There is no middle ground - you are either homosexual or heterosexual\nSexuality is a scale - it is possible to be somewhere in the middle",
       y = "") +
  theme_ggdist()+
  scale_y_continuous(
    labels = scales::percent_format(scale = 100))+
  theme(legend.position = "none",
        plot.title = element_text(face="bold"))+
  annotate("label", x = 1.2, y = .3, 
           label = "Men", size=4, color="#D60270", fontface="bold", hjust=0)+
  annotate("label", x = 1, y = .3, 
           label = "Women", size=4, color="#9b4f96", fontface="bold", hjust=0)+
  coord_flip()
ggsave("figures/appendix/FigureD3.png",dpi=500)

